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In the previous paper [Nenashev et al., arXiv:0912.3161 an analytical theory confirmed by nu- 
merical simulations has been developed for the field-dependent hopping diffusion coefficient D(F) 
in one-dimensional systems with Gaussian disorder. The main result of that paper is the linear, 
non-analytic field dependence of the diffusion coefficient at low electric fields. In the current paper, 
an analytical theory is developed for the field-dependent diffusion coefficient in three- and two- 
dimensional Gaussian disordered systems in the hopping transport regime. The theory predicts a 
smooth parabolic field dependence for the diffusion coefficient at low fields. The result is supported 
by Monte Carlo computer simulations. In spite of the smooth field dependences for the mobility and 
for the longitudinal diffusivity, the traditional Einstein form of the relation between these transport 
coefficients is shown to be violated even at very low electric fields. 

PACS numbers: 72.20.Ht, 72.20.Ee, 72.80.Ng, 72.80.Le, 02.70.Uu 



I. INTRODUCTION 

This paper represents a second part of our research 
dedicated to the theory of diffusion of hopping charge 
carriers biased by electric field in disordered systems with 
a Gaussian energy distribution of the energy of the local- 
ized states: 
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Here N is the spatial concentration of sites available 
for hopping transport and a is the energy scale of 
their density of states (DOS). The Gaussian DOS is as- 
sumed to apply for disordered organic materials, such 
as molecular ly doped and conjugated polymers and or- 
ganic glasses i£i2ii£i£i2 While in the previous paper— one- 
dimensional transport was considered, in the current pa- 
per we present results for two- and three-dimensional 
systems. This study has been on one hand stimulated 
by numerous experimental studies on organic disordered 
materials,— d ! 11 ! 12 ! 13 ! 14 ! 15 which claim the invalidity of 
the conventional form of the Einstein relation between 
the carrier mobility fj, and the diffusion coefficient D 



(2) 



On the other hand our study is stimulated by the lack 
of a concise theory for the diffusion biased by electric 
field in the hopping transport mode. It is, however, this 
transport mode that dominates the electrical conduction 
in disordered organic materials where transport is due 
to incoherent tunnelling of electrons and holes between 
localized states randomly distributed in space, with the 
DOS described by Eq. The transition rate 



between an occupied state i and an empty state j, sep- 
arated by the distance ry, is described by the Miller- 
Abrahams expression^ 
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where vq is the attempt-to-escape frequency. The energy 
difference between the sites is 
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where the electric field F is assumed to be directed along 
the Z-direction. The localization length of the charge 
carriers in the states contributing to the hopping trans- 
port is a. We assume the latter quantity to be inde- 
pendent of energy and we will neglect correlations be- 
tween the energies of the localized states, following the 
Gaussian-disorder-model of Bassler^^ 7 - 

The field-dependent diffusion in such systems in the 
3D case has so far been studied by computer simulations. 
Richert, Pautmeier, and Bassler performed Monte Carlo 
simulations in 3D and showed that the longitudinal diffu- 
sion coefficient D z is strongly dependent on the electric 
field and that the dependence is quadratic at such low 
fields that the mobility of charge carriers remains field- 
independent ,i2ii2i This is in contrast to the linear field 
dependence of the diffusion coefficient at low fields ob- 
tained by the exact analytical theory and by numerical 
calculations for ID systems in the preceding paper4 In 
order to clarify the nature of the field effect on the dif- 
fusion for 3D and 2D systems we suggest in the current 
paper an analytical theory for the field-dependent diffu- 
sion coefficient in the hopping regime for such systems. 
This theory confirms the conclusion of Richert, Paut- 
meier, and Bassle r 17 : 18 about the parabolic field depen- 
dence of D z at low fields. Furthermore, our theory gives 
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explicit analytical expressions for the combined effects 
of the electric field and temperature on the hopping dif- 
fusion coefficient. These expressions predict a violation 
of the traditional form of the relation between \i and D 
given by Eq. ([2]) at very low electric fields. Since the tem- 
perature effect on the field-dependent diffusion has been 
so far left out of the scope of computer simulations ; 17 ! 18 
we perform here a Monte Carlo study of the field- and 
temperature-dependent diffusion in the hopping regime 
for 2D and 3D systems in order to check the results of 
the analytical theory. Our computer simulations for the 
2D and 3D cases presented in Sec.|TT]support the analyti- 
cal theory. Furthermore, numerical results evidence that 
while the carrier mobility is stable with respect to differ- 
ent realizations of disorder, the diffusion coefficient ex- 
periences significant fluctuations from one realization to 
another, even for systems containing millions of localized 
states. The reasons for such different behavior between fj, 
and D is clarified and the estimates for the system size, 
which is necessary to obtain stable values of D, are given 
in Sec. Iff! CI The latter estimates show that previous nu- 
merical simulations in the literature were performed on 
rather small systems, insufficient for obtaining reliable 
results for the field-dependent diffusion coefficient at low 
temperatures. 

In developing the analytical theory for the field- 
dependent diffusion coefficient for the hopping transport 
regime we rely on the analogy between the hopping trans- 
port mode and the multiple-trapping mode i 19 ' 20 i 21 ' 22 In 
Sec. IIIII we present a general solution of the problem. 
In Sec. IIVI the analytical expressions from Sec IIIII are 
compared to simulation results. Concluding remarks are 
gathered in Sec. [V] 



II. MONTE CARLO SIMULATIONS 

In this section we study the effects of field and tem- 
perature on the diffusion coefficient by Monte Carlo sim- 
ulations in more detail than it has been done previously, 
aiming at a comparison with the analytical theory de- 
scribed in the following sections. Furthermore, we study 
in detail the role of the system size on the simulated re- 
sults and show that in order to get reliable results for 
the field-dependent diffusion at reasonably low temper- 
atures one needs to perform simulations on enormously 
large systems. 

The system is modelled as a lattice of L 2 or L 3 sites 
with lattice constant d and the site energies £j chosen 
randomly according to Eq. ([T]). Periodic boundary con- 
ditions are applied in all directions. Hops inside a square 
of 7 x 7 sites (2D case) or a cube of 7 x 7 x 7 sites (3D 
case) , centered at the starting site are allowed. The simu- 
lation proceeds as follows. A packet of n non-interacting 
carriers is allowed to move in the lattice until a fixed time 
t has passed. The mobility /i is calculated from the av- 
erage distance that the charge carriers have moved along 
the direction of the field, while the longitudinal diffusion 



coefficient D z is calculated from the width of the carrier 
packet: 



Further details of the simulation algorithm and our im- 
plementation of it are given in Appendix [X] Care was 
taken about the necessary size of the simulated system 
in order to avoid finite-size effects. The corresponding 
number of sites was 700 3 in the 3D case and I5000 2 in 
the 2D case. 

Simulation results for the diffusion coefficient both 
along and perpendicular to the electric field are shown 
in Fig. [TJ together with the mobility (scaled with kT/e). 
The localization length was a = Q.2d. A packet con- 
sisting of 1000 charge carriers was simulated for each 
data point, and the simulation results were averaged over 
five different realizations of disorder. At extremely low 
fields it is seen that all three of the plotted quantities are 
equal, which means that Einstein's relation, Eq. is 
valid. With rising magnitude of the electric field the lon- 
gitudinal diffusion coefficient increases drastically, while 
the mobility and transversal diffusion coefficient remain 
field-independent up to much higher fields. The solid line 
shows a fit to the longitudinal diffusion coefficient by a 
square trial function 

D Z (F,T) = A(T)F 2 + D {T). (6) 

It is seen that the square function well fits the data in 
agreement with the results of previous simulations^ 

In order to study the effect of temperature on the field- 
dependent diffusion coefficient, the simulations were re- 
peated for different temperatures. The results for the 
longitudinal diffusion coefficient in 3D are collected in 
Fig. [2] The data are fitted with the square trial functions 
©, shown in the figure by solid lines. The fit is good for 
all temperatures. The decrease of D Z {F) at very high 
fields for the two highest temperatures is due to the triv- 
ial saturation effect well-known for the one-dimensional 
random-energy models We focus in the following on the 
field dependence of D z at field magnitudes lower than the 
one at which D z starts to decrease with increasing field. 

The temperature dependences of the coefficient A in 
the quadratic field-dependent term for the longitudinal 
diffusion coefficient D z , that were obtained in the 2D 
and 3D cases by fitting the simulation results in Fig. [2] 
by Eq. ^ are shown in Fig. along with the field- 
independent term Dq. Analytical expressions for A will 
be derived in Sec. IIIII 

When performing computer simulations for the field- 
dependent diffusion one should be cautious with the 
choice of simulation parameters. The field-induced spa- 
tial spreading of the carrier packet is caused by trapping 
of some carriers onto localized states deep in energy, while 
the other carriers continue their motion in shallow states 
being biased by the electric field. In order to obtain re- 
liable results for the field-dependent diffusion coefficient, 
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FIG. 1: The diffusion coefficient along and perpendicular to 
the field, and the mobility fi scaled with kT/e as a function of 
the field strength, (a) in 3D, and (b) in 2D. At small fields all 
three quantities are equal, implying that the Einstein relation 
is valid. 
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FIG. 2: The longitudinal diffusion coefficient Dz for hopping 
in 3D, as a function of the applied electric field F. The solid 
lines show the best fit to the square function ((6]). 
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FIG. 3: The temperature dependence of Do, the diffusion 
coefficient at F = and A, obtained from Fig. [5] for the 
localization length a = 0.2d. 



one should guarantee the presence of such deep-in-energy 
states in the simulated system. Since the DOS given by 
Eq. dl]) rapidly decreases for the deep-in-energy states 
and hence the sites with deep energies are rare, one has 
to simulate large systems. The simulation time t must 
also be chosen so large that the charge carriers have time 
to visit the deep traps, which control the diffusion. In our 
simulations t was chosen for each temperature so that the 
typical number of hops for each carrier was at least 5 TO 7 . 

The importance of the deep-in-energy states for the 
field-dependent diffusion is demonstrated in Fig.[4j where 
the values of D z (at the field eFd = 0.5a) and of /i (at 
low fields) for the temperature kT — 0.33<7 are given 
when calculated with a cut-off of the DOS below some 
energy E c . For these calculations the normalization of 
the DOS was kept, while the states with energies below 
E c were excluded from the simulation. The mobility is 
almost unaffected by the cutting as long as E c < —4a, 
while the diffusion coefficient drastically decreases when 
sites with much smaller energies, around — 6a, are re- 
moved. The result for the mobility /x is not surprising. It 
has been predicted in the analytical theor y 21 ' 22 that the 
hopping mobility in the Gaussian DOS for the diluted set 
of carriers is determined by sites with energies in the close 
vicinity of the average carrier energy e av = —a 2 /kT. For 
kT = 0.33cr, this energy is e av ~ — 3cr, which explains the 
data for /i in Fig.|U The result for the diffusion coefficient 
in Fig.[4]shows however that rare sites with even lower en- 
ergies than e av cause the strong dependence of the diffu- 
sion coefficient on the electric field. In Sec. MI Cl we show 
that the most important sites for the field-dependent dif- 
fusion have energies around e* = —2a 2 /kT, which is 
much deeper than e av . This result also agrees with the 
data shown in Fig.H for kT = 0.33cr, e* w -6a. The dif- 
fusion coefficient changes drastically when the few sites 
with energies below — 5.5a are removed, while the mo- 
bility starts changing only when E c is in the vicinity of 
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FIG. 4: Diffusion coefficient D z (at the field eFd = 0.5a) 
and mobility /i (at low field) in a 3D system where sites with 
energies below E c are absent. The temperature is kT — 0.33<r. 



In our simulations a lattice of 700 3 sites has been used. 
There are typically no sites with energies below —6a in 
such a system. Therefore, the D Z (F) data for kT = 0.33cr 
cannot be considered as reliable in the whole simulated 
range of electric fields. However, for kT = 0.4c, the 
lowest temperature considered, e* = — ha. Energies 
around —5a are present in our lattice, and thus results 
for kT > 0.4cr can be considered as reliable for the lattice 
of 700 3 sites. It seems extremely lucky that in previous 
simulations with the lattice of just 70 3 site o 17 ' 18 mean- 
ingful results were claimed for kT = 0.33<r. 

The analytical results in the preceding paper— were 
obtained for nearest-neighbor hopping, while the simula- 
tions above allowed also longer hops. To exclude this dif- 
ference between the models as the cause of the different 
field dependencies, the simulations for two- and three- 
dimensional systems were repeated with only nearest- 
neighbor hopping allowed. The mobility and diffusion 
coefficient obtained in this case were somewhat lower 
than those for the variable range hopping, however the 
parabolic shape of the field dependence for the diffusion 
coefficient did not change. 

Our numerical results can be qualitatively summarized 
as follows: (i) the diffusion coefficient along the electric 
field depends parabolically on the field strength F at 
low fields; (ii) the diffusion coefficient perpendicular to 
the field is field-independent in the range of fields where 
Ohm's law is fulfilled; (iii) the field-dependent part of 
the diffusion coefficient rapidly decreases with increasing 
temperature; (iv) the field-dependent part of the diffu- 
sion coefficient is very sensitive to sites which energies 
are lower than the mean carrier energy. 

The discussion above, in particular the latter state- 
ment on the decisive role of rare sites with very deep 
energies that can hardly be found in the finite simula- 
tion arrays, raises the task of developing an analytical 



theory for the diffusion process in the hopping regime 
enhanced by an electric field. There is no such theory in 
the literature so far. The only relevant theory is the one 
developed by Rudenko and Arkhipov for band transport 
in materials with traps^ Although this theory predicts 
a parabolic field dependence of the longitudinal diffu- 
sion coefficient, it cannot be directly applied to hopping 
transport, because it operates with quantities that are 
specific for band conductivity (for example, the effective 
density of states in the conduction band). Therefore it 
is necessary to develop a theory for hopping transport, 
particularly because Monte Carlo simulations suffer from 
finite-size effects as described above. We will give such a 
theory in Sec. IIIII 



III. ANALYTICAL RESULTS 

A. Approximation of independent jumps 
(general consideration) 

The origin of the field-dependent diffusion can be un- 
derstood qualitatively with the aid of a spatio-temporal 
picture of the carrier distribution sketched in Fig. [5] The 
small dots indicate the scatter of carriers after some def- 
inite number of jumps, assuming that all carriers start at 
the same time t = from the same point. To get the spa- 
tial distribution of carriers at some definite time t* , one 
may "project" these dots from the starting point to the 
line t — t*. The direction of "projecting" is determined 
by the drift velocity. When the drift velocity is equal to 
zero (no electric field, the left part of Fig. [5]), the spatial 
distribution of carriers at t = t* docs not depend on the 
scatter of times spent by the carriers in order to perform 
a fixed number of jumps. On the contrary, in the case of 
drift caused by an electric field (the right part of Fig. , 
this scatter in times "projects" into the line t = t* , which 
gives rise to the broadening of the spatial distribution at 
t = t*. This broadening is the reason for the enhance- 
ment of the diffusion coefficient due to electric field. 

The above consideration shows that fluctuations of 
the durations of jumps are responsible for the field- 
induced diffusion. These fluctuations can be especially 
pronounced for systems with a broad distribution of site 
energies, because jumps from energetically deep sites to 
transport sites demand exponentially long times. 

Let us start from a general form of the quantitative 
description for the field dependence of the diffusion co- 
efficient. Our consideration is based on the assumption 
that successive jumps are statistically independent, i. e. 
the diffusion process is Markovian. The latter means 
that increments in carrier coordinates at a given jump, 
as well as the time interval from the preceding jump till 
the given one, do not depend on the carrier prehistory. 
This assumption definitely does not hold for the one- 
dimensional hopping transport considered in the previ- 
ous paper £ In the ID case, the probability of returning 
to an already visited trap is not negligible, and hence 
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FIG. 5: A sketch of spatio-temporal distribution of carriers 
after some definite number of jumps. 



the consequent jumps must be correlated. Therefore our 
analytical consideration based on the assumption of the 
statistical independence of the successive jumps present 
below can be valid only for 2D and 3D cases since it is 
reasonable to assume that in the latter cases the carrier 
trajectories are non-returning. 

Let t* be some fixed moment of time (t* is large com- 
pared to the mean time of a jump); x* , y* and z* be the 
carrier's displacements along the axes X, Y and Z during 
the time interval [0, t*]; x\, y±, Z\ and t\ be coordinate 
displacements and a time increment related to only one 
jump, respectively; x n , y n , z n be the displacements after 
n successive jumps. Angle brackets will denote averaging 
over jumps performed by different carriers, equivalent to 
the averaging over successive jumps of one carrier. 

We start with the simple case of zero electric field. 
Since there is no drift, the expectation values of x*, y* 
and z* vanish, and one obtains for the diffusion coeffi- 
cients D x , D y , D z 



D x = lim 

t* — >oc 



2t* 



-, Dy = lim 



2t* 



D z = lim 



t*-»°o 2t* 



(7) 

For large t* , the displacements x* , y* , z* are approx- 
imately equal to xn, Vn, z n ~ displacements after TV 
jumps, where A^ = t* / (h) is the mean number of jumps 
during the time t*. Consequently, 
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In the limit t* — > oo this approximate equality becomes 
exact, and one gets 
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Let us now consider the case of a finite electric field 
along the Z-axis. Since the expectation values of x* and 
y* are still zero, all the considerations above remain valid, 



and Eq. ([9]) remains correct with respect to D x and D y . 
One cannot obtain D z by literally the same way, because 
(z*) ^ 0. Instead, one may apply this argumentation to 
a variable z = z — vt, where v = (z\) / (h) is the drift 
velocity. Since (z*) = (zn) = (zi) = 0, one gets 



D z = D z = lim 



(z~* 2 )_ (zl) (zl)-2v( Zl h)+vHtl) 



f^oc 2t* 



2(h) 



2(h) 



and finally 
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(11) 



Since the mean values (zi), (z^) and (zih) depend on 
the electric field, Eq. (jTTJ) describes the field-dependent 
diffusion along the field direction. On the other hand, 
there is no reason for the mean values (xf ), (yf) and (h) 
to be field-dependent in small electric fields. Therefore, 
according to Eq. ©, the transversal diffusion coefficients 
D x and D y are expected to be constant inside the Ohmic 
regime. 

Let us discuss the shape of the dependence D z (F) near 
the point F = 0. Since (z\) is proportional to the electric 
field ((z\) = fiF(t\), where jU is mobility), the third term 
in the r.-h. side of Eq. (|11[) gives a contribution to D Z (F) 
that is quadratic in F. In the first term, 



(z\) = ( Zl ) 2 + a 2 [ Zl ] 



(12) 



where a[a] denotes the standard deviation of a random 
variable a. Since cr[z\\ is not sensitive to the electric field, 
the first term in the r.-h. side of Eq. pip is a sum of a 
constant and a term quadratic in F. The behavior of the 
second term in the r.-h. side of Eq. (JTTJ) depends on the 
symmetry of the system. If in the absence of electric field 
the directions Z and —Z are equivalent, then (zih) = 
at F = 0. In this case one expects that (zih) ~ F and 
consequently the second term in the r.-h. side of Eq. (jTTJ) 
is quadratic in F. However if the positive and negative 
directions along the Z-axis are non- equivalent, (ziti) can 
be non-zero at F = 0, which gives a contribution to the 
diffusion coefficient linear in F. 

In summary, if the directions Z and —Z are equiva- 
lent, which we will assume in the following, the longitu- 
dinal diffusion coefficient D z at small fields is described 
by Eq. © where Do and A are field-independent coeffi- 
cients. 



B. Multiple-trapping band conductivity 

Let us first apply our approach to the multiple- 
trapping (MT) model of band conductivity. The model 
includes processes of capture of free carriers by traps, 
emission of trapped carriers, and free motion (Brownian 
motion plus drift in the electric field) of carriers in the 
band (Fig.[S]). The band is characterized by the effective 
density of states N c , the mobility /// and the diffusion 
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FIG. 6: Multiple-trapping conductivity. 



coefficient Df of free carriers. The traps are character- 
ized by the density of states g(e) and the capture rate 
(for unit free carrier concentration) c(e). The energy e is 
counted from the band edge. 

The problem of the field-effect on the diffusion coeffi- 
cient in the MT model has been considered by Rudenko 
and Arkhipov{2£ who treated the evolution in time of the 
one-dimensional carrier distribution function. In Ref. [23l . 
transport was assumed to take place in quasi-equilibrium, 
i. e., the distribution function was assumed to change 
slowly in comparison with the rate of exchange between 
traps and the conduction band. Our approach is free of 
the assumption of quasi-equilibrium transport. Besides, 
our consideration provides information not only about 
the longitudinal diffusion coefficient, but also about the 
transversal one. 

Instead of examining the carrier distribution function, 
we will follow the motion of a single carrier and consider 
the statistics of this motion. In order to use the expres- 
sions © and we represent the motion of a carrier 
as a sequence of "jumps" , each of them beginning at the 
moment of escape from a trap. Successive "jumps" are 
statistically independent, because the processes of cap- 
ture take place in different points in space and there is 
no correlation between them. 

Each "jump" of a carrier consists of two contributions: 
free motion and sitting on a trap. Let t\f be the time 
of free motion (between emission and capture), and tit 
be the time of being trapped (between capture and emis- 
sion); t±f and tit are independent random variables. We 
denote their expectation values (ti/) and (tit) as Xi and 
T2, respectively. Since Xi is the lifetime of free carriers, 
its reciprocal value is the capture rate: 



xr 



J c(e)g(e)de. 



(13) 



The variable tit obeys an exponential distribution. 
Hence its mean square (t\f) is equal to 2X 2 . 

Collecting all this information, one gets the following 
statistical expressions for a single "jump" (electric field 



F is directed along Z): 
(h) =Ti+T 2 , 

(if) = (t\ f ) + 2{t lf t u ) + (ti) =2Xi(Xi + X 2 ) + (tit), 
(xi) = (yi) = 0, (14) 
(21) = HfF{t X f) = HfFT x , 
(xi) = (y 2 i)=2D f (tx f )=2D f Tx, 
(4) = ^}F 2 (t\ f ) + 2D f (tx f ) = 2 M 2 F 2 X 2 + 2D f Ti, 
(zitx) = (zxtxf) + (zxtxt) = 2v f FT? + v f FTiT 2 . 
Using these mean values, one can evaluate the mobility: 



/' : 



Xi 



F(ti) ^ J Ti+T 2 
the transversal diffusion coefficient via Eq. 

Xi 



(15) 



D x = D„ 



(*?) 



D 



2(h) / X 1 +X 2 ' 
and the longitudinal diffusion coefficient via Eq. (fTTj) 



(16) 



D 7 =D 



Xi + Xj 



2(T 1 + T 2 ) a 



(17) 



The parameters Xl, T 2 , and (i 2 t ) are governed by the 
capture cross-sections and emission rates of the traps. 
For moderate electric fields, these cross-sections and rates 
can be regarded as field- independent, and consequently 
one can use equilibrium values for Xi, T 2 , and (t{ t ). This 
gives an opportunity to simplify Eqs. (|15p - (|17[) . The ra- 
tio a = Xi/(Xl + X 2 ) contributing to these equations 
is simply the fraction of free carriers in the equilibrium 
state. The distribution of the dwell times r of a carrier at 
some individual trap is exponential with the mean value 
(r) = 1/r-f, where Tj is the emission rate. Consequently, 
the mean square of this dwell time is (r 2 ) = 2TJ 2 . The 
mean value (i 2 t ) is a weighted average of values (r 2 ), 
where the probability of visiting a trap serves as the 
weight. Therefore, 



(tit) 



2Tf(e) P (e)de, 



(18) 



where p(e)de is the probability of visiting (at a given 
"jump") a trap with an energy in the range [e; e + de]: 

^'jM^ sT ' c(£)s{s)d£ - (i9> 

Expressing the escape rates through capture cross- 
sections: 



r T (e) = c{s)N c e' 



e/kT 



one obtains (tl t ) = 2TiI, where 



I = 



c(e)N c 



c(e)g(e)ds. 



(20) 



(21) 
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Finally, Eqs. (fT5l) - (fT7|) get the following form: 

H = afi f , D xy = aD fl D z = aD f + a 3 In 2 f F 2 , (22) 

where a is the equilibrium fraction of free carriers, and 
I is the integral defined by Eq. pTjl . 

The expressions for fi and D z are the same as the ones 
obtained by Rudenko and Arkhipov. 23 However, our con- 
sideration provides some new information about diffu- 
sion in the multiple-trapping conductivity regime. First, 
Eq. (|22| shows that the transversal diffusion coefficient 
is field-independent. Second, we have proven that the 
equations for /i and D z are exact (provided that emis- 
sion/capture probabilities are not influenced by electric 
field) , whereas in Ref . [23j they were obtained only under 
so-called "quasi-equilibrium conditions" . 



durations of these contributions; Ti and T2 denoting the 
mean values (t\f) and (tit), respectively (where averaging 
is over successive "jumps"); [if and Df be the mobility 
and the diffusion coefficient of "free" motion (energeti- 
cally far above the traps). With these notations one can 
follow the same derivation as in the previous subsection 
and see that Eqs. (fT5|) -(fT7 |) are still valid also for hopping 
transport, ft is convenient to rewrite these equations: 

M = Mo, D x = D y =D , D Z = D + AF 2 , (23) 

where 

MO = P>f rp , rp ) A) = Uf — — , 



C. Hopping transport 

Now we will apply Eq. (fTTj) to two- and three- 
dimensional hopping transport in a system with a Gaus- 
sian DOS, in a manner very similar to our consideration 
of the multiple-trapping conductivity. 

Our analysis is based on the following hypothesis: the 
field dependence of the diffusion coefficient in 2D and 
3D hopping is mainly due to very rare and energetically 
deep sites. We will call them "traps". The traps are 
rare in two senses: first, the typical distance between the 
traps is large in comparison to the inhomogeneities of 
the mobility; second, the probability of being trapped is 
small, so that the traps do not affect the carrier mobility. 
"Energetically deep" means that the energies of the traps 
are far below the mean energy of the carriers. We will 
show below that this hypothesis provides a reasonable 
description of the Monte Carlo simulation results on the 
field-dependent diffusion. 

As in the previous subsection, we consider the mo- 
tion of a carrier as a sequence of "jumps" , each begin- 
ning when the carrier enters a trap, and ending when 
the carrier enters another trap. The durations of dif- 
ferent "jumps" are not correlated. The same is true for 
the carrier displacements. Indeed, the three-dimensional 
Brownian motion is non-returning. This property guar- 
antees that the carrier always visits a new trap, and that 
the trajectories of its motion between different traps do 
not overlap. Hence there are no reasons for correlations 
between successive "jumps" . This is the point where di- 
mensionality is important. In one dimension, the proba- 
bility of visiting a previously visited trap is not negligible. 
Consequently the "jumps" must be correlated. Since the 
"jumps" in the 3D case are not correlated, one can obtain 
the mobility \i and the diffusion coefficient D from the 
statistics of "jumps" using the method of Sec. MI Al For 
the 2D case we also use the assumption of independent 
"jumps" . 

As in the case of multiple-trapping conductivity, each 
"jump" consists of two contributions: "free" motion of a 
carrier and "sitting" on a trap. Let t\f and t\ t be the 



A : 



Mo 



2 (tl) 



2(Ti + T 2 ) ' 



(24) 



In the regime of ohmic conductivity, the values of /if, 
Df, Ti, T 2 and {t\ t ) can be regarded as field- independent 
because a sufficiently small electric field does not signif- 
icantly perturb the probabilities of capture and release. 
Therefore one can neglect a possible dependence of /xo, 
Do and A on the electric field. In the following we will 
use the zero-field values for these quantities. 

Let us calculate the coefficient A. For convenience, we 
will consider the system as a large but finite one (with 
periodic boundary conditions, to allow drift). Then one 
can obtain the following expressions for Ti, T2 and (t\ t ): 



Ti = 



T s 



,PtT. 



CSC,t 



pt ■ 



(t 2 lt ) = 2T 1 J2PtT. 



1 

esc.t > 



(25) 



(26) 



(27) 



where the index t runs over all traps, pt is the probability 
that a carrier is at site t, and r esCjt is the rate of escaping 
from the trap t. 

Eq. (|25|) results from a consideration of the carrier flow 
from/to traps. Namely, the flow of carriers out of traps 
is equal to ^2 tP t^esc,t] the flow into traps is Tf 1 . In 
equilibrium, these flows are equal to each other, which 
gives Eq. ([2"5]) . 

In order to obtain Eqs. (|26[) . (|2T|) . we introduce the 
probability Vt that trap t will be the next visited trap. 
Let us consider the balance of flows from/to trap t. Flow 
from this trap is equal to ptr esc ,t- Flow to it is VtT^f . 
Therefore 



V t = T lPt T c 



(28) 
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The mean time of being captured at trap t is r~ c ( . To 
obtain T 2 , we average these mean times with correspond- 
ing weights Vt- 



Tl = reseat . 
i 



(29) 



Substituting here Eq. ([25]). one obtains Eq. (j2"6")l . Analo- 
gously, the mean square of the time of being trapped at 
site t is 2r~ 2 t (the factor of 2 comes from the exponential 
distribution of dwell times). Again, we take a weighted 
average to obtain (t\ t ): 



(30) 



that, together with Eq. (28), provides Eq. (gTJ). 

The substitution of Eqs. [p5 |) -(|2? f into Eq. ([21]) gives 
the following expression for A: 



(31) 



Let us try to simplify this expression. The sum in the 
denominator is the probability that a carrier is at some 
trap. Due to the rarity of traps, this probability is small, 
and it can be neglected. The summation in the numera- 
tor will be taken over all sites: 



2 

Mo 



(32) 



We will use curly braces to denote the summation in 
which index s runs over all sites: 



P s a s , 



(33) 



where a is any quantity specific for sites. It is obvious 
that curly braces mean averaging over an ensemble of 
particles (or time averaging, which is the same for finite 
system and long enough time). Thus, 



A w fig {r es *} . 



(34) 



Then, the rate of escape r esCjS is related to the sum 
^2 s '^s Fas' of transition rates from site s to any other 
sites: 



E^ 



(35) 



where n eS c.s is the mean number of "attempts" to es- 
cape from site s (including the successful one). These 
"attempts" are events of carrier's hopping out of site s 
before the carrier moves so far from this site that it com- 
pletely forgets the prehistory related to this site. Let us 
denote the sum ^] s ,^ s T ss i as tj . Then one can rewrite 
A as 



A ~ /ig \n esc i\ 



(36) 



or, introducing the averaged number of escape attempts 



^csc i 



A ~ n csc Hq {t} . 



(37) 



Calculating {t} implies averaging over site energies (s s ) 
and over quantities related to the neighborhood (energies 
of neighbors e s > and distances to them r ss /). It is possible 
to treat these two kinds of averaging separately, using 
our assumption that "optimal" traps are deep in energy. 
Since the traps are deep, each hop from a trap is upward 
in energy. Therefore the Miller- Abrahams rates, Eq. , 
for such hops can be written as 



r. 



vq exp 



2r, 



kT 



(38) 



The dependence of this expression on e s has the form 
of a factor exp(e s /fcT). Hence, the quantity t s = 
Tss') 1 can b e factorized (only for deep sites s) as 



t s = t s cxp(-e s /A;T), 
where t s does not depend on e s : 



E ex p 



2r, 



ffl 
kT 



(39) 



(40) 



Since there are no correlations between exp(— e s /kT) and 
one can average these quantities separately: 



{t} » { t 



-es/kT T x _ 



-Es/kT 



} 



(41) 



Here r is the mean value of r s (an arithmetical average 
over all s). Now, the calculation of {e~ 6s ^ kT } is simple. 
Remember that p s is an equilibrium probability of finding 
a carrier at site s: 



-es/kT 



, /kT ' 



Therefore 



-e s /kT\ — 



-e B /kT _ 



-2e s /kT 



^ 6 _ Je- 2 "/ kT g(e)de 

£ e -e s ,/fcT " f e -c/ kT g{e)de ' 



(42) 



(43) 



(44) 



For a Gaussian DOS g(e) given by Eq. |T]), one obtains 

3a 2 



{e -es/kT } = exp 



(- 



Hence 



{t} w t exp 



\2(kT) 2 

f 3a 2 
WkT)" 



(45) 



(46) 
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Substituting this result into Eq. (|37|) . one obtains 



.4 



3a 2 



(47) 



\2(kT) 2 i 

This equation contains r. The proper way to calculate 
r is the numerical one: to generate at random a large 
enough number of neighborhoods of the given site (i. e. 
sets of energies e s > of neighboring sites and distances r ss > 
from the given site s), to calculate r s for each neigh- 
borhood according to Eq. (|4H)) . and then to average the 
results for r s . In the next section we compare Eq. ([47]) 
to the results of Monte Carlo simulations from Sec. [TTJ 

Let us now reveal which traps are "optimal" in the 
sense that they determine the diffusion coefficient. The 
integrand in the numerator of Eq. (|44[) has a sharp peak 
of width a at e = —2a 2 /kT. Hence the "optimal" traps 
are sites with energies in the range —2a 2 /kT ± a. 

It is now possible to justify our assumption about 
the possibility to distinguish between the transport sites 
and the traps. We restrict ourselves to the case of 
low temperatures, kT <C a. The typical energy of the 
"optimal" traps, —2a 2 /kT, is significantly lower that 
the mean energy of carriers {e} = —a 2 /kT (the differ- 
ence is much larger than kT). Herewith the assump- 
tion that traps are deep in energy is fulfilled. The 
probability for a carrier to have an energy lower than 
-2a 2 /kT is - exp(-er 2 /2(fcT) 2 ). This value is much 
smaller than unity. Therefore neglecting the sum in the 
denominator of Eq. (|31[) is justified. Let us now find 
the typical distance L t between "optimal" traps. Their 
concentration can be estimated as the concentration of 
sites with energies lower than —2a 2 /kT , which is about 
~ N exp(— 2a 2 /(kT) 2 ), where N is the total concentra- 
tion of sites. Hence L t ~ 7V~ 1/3 exp(2cr 2 /3(/cT) 2 ). We 
should compare it to the characteristic size L M of in- 
homogeneities of the mobility. One can estimate this 
size as a typical distance between sites that are im- 
portant for the mobility. These are sites with ener- 
gies close to the mean energy of carriers ; 21 ' 22 —a 2 /kT. 
Their concentration is about N exp(— a 2 /2(kT) 2 ). Hence 
~ 7V~ 1 / 3 exp(cr 2 /6(fcT) 2 ). Herewith we obtain the 
strong inequality L t ^> L^, which implies that it is safe 
to describe the motion of the carriers between traps by 
means of a constant mobility fj,f and a diffusion coeffi- 
cient Df. 

All of the above considerations were based on the 
Boltzmann statistics for the charge carriers. Let us dis- 
cuss the effect of the carrier concentration n on the 
obtained results. The usage of Boltzmann statistics is 
justified if the Fermi level Ef is far below the ener- 
gies of the sites that make a major contribution to the 
mean value {t}. As described above, the main contribu- 
tion to {t} comes from sites with energies in the vicin- 
ity of ~2(T 2 /kT. A simple calculation shows that car- 
rier concentration ridis corresponding to a Fermi level 
c F = -2a 2 /kT is 



where N is the concentration of sites. The theory con- 
sidered above presumes that n -C ^diff- In this limit, 
the diffusion coefficient does not depend on the carrier 
concentration. For larger concentrations, n 3> ridiff, sites 
with energies w —2a 2 /kT are essentially occupied. In 
the latter case they cannot efficiently capture the mov- 
ing charge carriers and therefore the contribution of such 
sites to the diffusion process (in particular to {t}) is sup- 
pressed. According to Eq. (f37|) . the field-induced part of 
the diffusion coefficient D Z (F) — D is proportional to 
{t}. Therefore, in the case n ^> UdiS, D Z (F) — Dq de- 
creases with increasing n. Detailed consideration of this 
concentration-dependent diffusion is beyond the scope of 
the present paper. 

We should also note that if the sample contains less 
than ~ exp(2c 2 /(fcT) 2 ) sites, then the average number 
of "optimal" traps in the sample is less than 1. One 
should therefore expect large sample-to-sample variations 
of D Z (F) — D in this case. 



IV. COMPARISON OF MONTE CARLO AND 
ANALYTICAL RESULTS 

In this section we discuss relations between the simu- 
lation results described in Sec. |TT] and the analytical ex- 
pression (|4"7|) . 

There are three quantities in Eq. (|47[) that are not 
input parameters of the model: the zero-field mobility 
/io, the averaged number of escape attempts n^, and 
the value of r. For the mobility, we will use the values 
taken from our simulations. There is no need for a special 
discussion of the mobility, because it is a well-studied 
property of the transport model considered here i 1 ' 21 ' 22 

The value of r was obtained numerically, as explained 
in the previous section with respect to Eq. (I47[) . The 
calculation of r is relatively inexpensive (as compared to 
Monte Carlo simulations of the diffusion) and it can be 
performed over a wide range of model parameters. For 
the lattice model used in Sec. HH we have found that the 
simulated dependence of r on model parameters is well 
fitted by the following expressions in two (tzd) and three 
(t~3d) dimensions: 



r\2D 



t 3D = cxp (\ 3D {j^) 1 + V3d \ , 



(49) 



(50) 



where the coefficients A and 77 depend on the ratio of the 
localization radius a to the lattice parameter d: 



X 2D = -1.12 (a/d)°- 2 + 0.76, 
V 2D -2.55 {d/a)°- 9 - 2.4, 
A 3 d = 0.27 log(d/o) - 0.74, 
ria D = 2.92 (d/a)°- 85 - 3.28. 



(51) 
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The error of fitting does not exceed 2.5% within the range 
of parameters 0.2 < kT/a < 1 and 0.2 < a/d < 0.5. 

Extensive computer simulations are necessary to de- 
termine the exact value of n osc . Without simulations one 
can only claim that n csc is a number larger than unity, 
though not exponentially large as a function of model 
parameters kT/a and a/d. Indeed, it follows from the 
concept of transport energy , 20 ' 21 ! 22 that the energy of a 
carrier just after its hop from a trap is close to the trans- 
port energy. It means that there is no energetic barrier 
for moving further away from a trap after the first hop. 
Therefore it is natural to suppose that, when a carrier 
has jumped out of a trap the probability of escaping i. e. 
(n osc ) _1 is comparable to 1. 

In order to learn more about n esc , we used the values 
of /io and A obtained by Monte Carlo simulations. Sub- 
stituting the data shown in Fig. [3] into Eq. (j4"Tj) , we have 
extracted n osc for different temperatures in the range 
0.33fi < kT < 0.7a. The localization radius a was chosen 
equal to 0.2<i. We have obtained that varies in the 
range from 4.5 to 8 in the 2D system, and in the range 
from 1.6 to 3.1 in the 3D system. Simulations have also 
shown that the values of do not significantly change 
when the localization radius is changed from 0.2<i to 0.5d. 
One can see that values of n csc found by simulations are 
indeed reasonable: they are larger than unity but of that 
order. 

Therefore, the analytical consideration presented in 
Sec. IIII CI reduces the problem of predicting the field ef- 
fect on diffusion to (i) the evaluation of the zero-field 
mobility /in, and (ii) the evaluation of a coefficient ri osc , 
which is a slowly varying function of the system param- 
eters kT/a and a/d. The mobility can be obtained ei- 
ther from computer simulations (which is a much easier 
task than simulations of the diffusion), or from the ana- 
lytical theorji 2 ^ and experiments. For the sake of self- 
consistency, we rely here on the simulation data for the 
mobility. With respect to n esc , the comparison between 
numerical results and the theory shows that one can con- 
sider this coefficient as a constant and nevertheless ob- 
tain results that agree with simulations. The appropriate 
choice of this constant for a = 0.2 d is 

— -I 6 ' (2D) ' (52) 
Ucsc ~ \ 2.3 (3D). 



In Fig. [Jj Monte Carlo results for the coefficient A, 
which describes the field-induced diffusion according to 
Eq. ([B]), are compared to the analytical expression (|4T)l . 
The coefficient n^: in Eq. (|4"T)) is set to a constant accord- 
ing to Eq. (|52p . It is evident that in the framework of the 
simplifying assumption of constant n csci the analytical 
theory correctly reproduces the shape of the temperature 
dependence for the field- induced diffusion coefficient. 



1 


1 1 


MC simulation, 2D O 


o 




MC simulation, 3D • 








# 
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FIG. 7: Comparison of Monte Carlo results for coefficient 
A (symbols) with Eq. (|47p (lines) for 2D and 3D transport. 
Localization length a — 0.2d. Values of r and taken 
from Eqs. fl9) , ||5D) , ||52) . For no Monte Carlo results are used. 



V. CONCLUSIONS 

The main result of this paper is the development of 
an analytical theory for the field-induced diffusion in the 
hopping transport mode in 3D and 2D systems with the 
Gaussian DOS given by Eq. (fT]). At low electric fields, the 
field dependence of the longitudinal diffusion coefficient 
D Z (F) is parabolic as expressed in Eq. ([6]). The analyti- 
cal expression (|4"T]) gives the temperature dependence of 
the field-induced diffusion. Accompanying Monte Carlo 
simulations confirm the analytical results and show that 
the shape of the field dependence is parabolic. Together 
with the exact results of the previous paper— stating the 
non-analytic linear field dependence of the diffusion co- 
efficient in the ID case, our result resolves the long- 
standing puzzle for the reason of the discrepancies be- 
tween the analytic and non-analytic behavior of D Z (F) 
claimed in different studie a 18 ' 24 : It is the space dimen- 
sionality that is responsible for non-analytic (ID) and 
analytic (2D and 3D) dependences in D Z (F). 

Furthermore, our theory shows that the main contri- 
bution to the field-induced diffusion process comes from 
localized states with energies in the vicinity of —2a 2 /kT. 
The DOS parameter a in organic semiconductors is of 
the order of 0.1 eVi^i Therefore at room tempera- 
tures this energy —2a 2 /kT, which is decisive for D Z (F), 
is situated very deep in the tail of the DOS, around — 8a. 
This fact raises very severe demands to the size of the 
system in computer simulations, which aim at studying 
the field-induced diffusion in organic semiconductors. In 
order to have the decisive traps in a simulation at room 
temperature one needs approximately 10 16 sites. There- 
fore, such simulations cannot be considered suitable for 
studying the field-induced diffusion. For instance, the 
system size of 70 3 used in the previous simulation o 17 ' 18 
is suitable only at kT > 0.55cx. In contrast, our Simula- 
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tions carried out on systems with 700 3 sites give reliable 
results at kT > OAa confirming the developed analytical 
theory in a wide range of parameters. 

In all models the influence of the electric field on 
the mobility and on the diffusion coefficient increases 
with decreasing temperature. For the mobility this phe- 
nomenon has been accounted in the frame of the concept 
of the effective temperature ! 25 ! 26 The results of this pa- 
per show that for the longitudinal diffusion coefficient the 
concept of the deep traps should be used instead of the 
effective temperature. 
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Calculating the hopping rates ([3]) is very time consum- 
ing since the exponential function is expensive to com- 
pute. If we instead of the site energies £j store their 
exponentials, 



exp i 



kT) 



(A5) 



the hopping rates Tij can be computed more efficiently. 
We use the fact that only a small number of discrete dis- 
placements are possible for hopping in a lattice, when we 
restrict the length of the hops. Therefore the geomet- 
ric part of the hopping rate and the energy contribution 
from the electric field can be calculated and stored once 
for each displacement. In our case the hops are restricted 
to a cube of 7 x 7 x 7 sites centered at the starting site. 
For each displacement Ar, define the quantities 



PAr 



exp 



Ar 
-2 — 

a 



(A6) 



and 



APPENDIX A: MONTE CARLO ALGORITHM 

For a numerical study of the diffusion in hopping trans- 
port, an algorithm is needed that can efficiently simulate 
transport in large systems. Since the number of sites that 
can be treated in the simulation is limited by the available 
memory, we have chosen to study hopping transport on 
a lattice instead of a system with randomly placed sites. 

When the charge carrier is located at site i, the prob- 
ability that the next jump takes it to the site j is given 
by 



-L i 



(Al) 



where is the total rate of hopping away from site i: 



(A2) 



The time r that the charge carrier spends on the site i 
before hopping, (the "dwell time") is calculated as 



T/Tt 



(A3) 



where T for each hop is randomly generated with an ex- 
ponential distribution with unit variance. Which jump 
to perform is decided by picking a random number x be- 
tween and 1 from a uniform distribution, and finding j 
such that 



(A4) 



fc=i 



fc=i 



This ensures that each site k is selected with probability 
Pk ■ So far this is the standard Monte Carlo algorithm for 
hopping transport] 5 ' 17 ! 18 Below an efficient implementa- 
tion of this algorithm will be described. 



fAr 



exp 



zFAz 
kT 



(A7) 



where Az is the z-component of Ar. The hopping rate 
from site i to site j located at the position Ar relative to 
i can now be evaluated using 



y — z^o/OAr-min ( 1, — fAr 



(A8) 



One could also consider the storing of all hopping rates 
Tij , but since this greatly increases the amount of mem- 
ory needed by the simulation, it would restrict the size 
of the systems that can be simulated. We have found 
that a good balance between the simulation speed and 
the memory requirements is achieved by storing only the 
rate Ti and the quantity for each site, and to calculate 
the rates by Eq. I|A8[) during the simulation each time 
they are needed. 

One further optimization is to consider the jumps in 
order of increasing lengths, since shorter jumps typically 
have higher probabilities than longer ones. This ordering 
greatly reduces the number of hopping rates that have to 
be evaluated before the destination site j that satisfies 
Eq. (|A4|) is found. 

Each carrier was initially placed on a randomly chosen 
site in the lattice and then allowed to equilibrate before 
the measurements of time and displacements started. 
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